rm(list = ls())
library(ggplot2)
library(patchwork)
num <- 100
gm <- 0.95
eps <- 0.01

# eps冲击下技术水平的动态演进
# eps第0期为0，第1期为0.01，其他期为0
picdata <- data.frame(t = 1:num, lambda = NA)
picdata$lambda[1] <- 0
picdata$lambda[2] <- eps
for (i in 3:nrow(picdata)) {
  picdata$lambda[i] <- gm * picdata$lambda[i-1]
}
p1 <- ggplot(data = picdata, aes(x = t, y = lambda)) + geom_line() + theme_bw()
# ggsave('../lambda.png')

# eps冲击下资本存量、消费、产出、劳动时间、利率的动态演进
picdata$r <- picdata$H <- picdata$Y <- picdata$cons <- picdata$K <- NA
picdata$r[1] <- picdata$H[1] <- picdata$Y[1] <- picdata$cons[1] <- picdata$K[1] <- 0
for (i in 2:nrow(picdata)) {
  picdata$K[i] <- 0.9537*picdata$K[i-1] + 0.1132*picdata$lambda[i-1]
  picdata$cons[i] <- 0.5691*picdata$K[i-1] + 0.392*picdata$lambda[i-1]
  picdata$Y[i] <- 0.2045*picdata$K[i-1] + 1.4523*picdata$lambda[i-1]
  picdata$H[i] <- -0.2430*picdata$K[i-1] + 0.7067*picdata$lambda[i-1]
  picdata$r[i] <- -0.7955*picdata$K[i-1] + 1.4523*picdata$lambda[i-1]
}
p2 <- ggplot(data = picdata, aes(x = t, y = K)) + geom_line() +
  geom_line(aes(y = cons), linetype = 2) +
  geom_line(aes(y = Y), linetype = 3) +
  geom_line(aes(y = H), linetype = 4) +
  geom_line(aes(y = r), linetype = 5) + labs(y = '') +
  annotate('text',x = c(25,20,12,25,25), y = c(0.009,0.007,0.012,0,-0.0028), label = c('K','C','Y','H','r')) +
  theme_bw()
p1/p2
# ggsave('../irf.png')
